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We present an extended discussion of a recently proposed theoretical approach for off-resonance 
tunneling transport. The proofs and the arguments are explained at length and simple analogies and 
illustrations are used where possible. The result is an analytic formula for the asymptotic tunneling 
conductance which involves the overlap of three well defined physical quantities. We argue that 
the formula can be used to gain fresh insight into the tunneling transport characteristics of various 
systems. The formalism is applied here to molecular devices consisting of planar phenyl chains 
connected to gold electrodes via amine linkers. 

PACS numbers: 



Exciting new developments in molecular transport 
have lead to accurate single molecule measurements. A 
large part of the experimental data has be en confirm ed 
by several independent experimental groups. ^^ * '^ * ^ * ^*^ ! We 
are particularly interested in the data for molecules made 
of several repeating units or monomers, such as alkyl, 
phenyl, acene chains of various lengths. For such systems, 
we develop a semi-analytic theory of tunneling transport. 

The signature of the tunneling transport is exponential 
dependence, G=Gce~^^ , of the conductance G on the 
number of monomers N. In the past, the off-resonant 
tunneling transport was described and understood in 
terms of effective electrons tunneling through square 
barriers.l^ Such treatment works well as long as the effec- 
tive mass approximation remains valid at the Fermi level. 
However, many systems, in particular organic chains, dis- 
play large insulating gaps and flat bands and very often 
the effective mass approximation for these systems fails 
when one moves away from the band edges. The modern 
theory of tunneling transportPsESEl connects the tun- 
neling exponent /3 to the complex band structure of the 
chains, an approach that goes well beyond the effective 
mass treatments. In a past public at ion ,^^2! we have con- 
tributed to the picture by deriving an expression for Gc, 
the contact conductance. 

In this work, we extend our previous discussion of 
the off-resonant tunneling transport in periodic insu- 
lating chains. We carefully review our previous argu- 
ments, extend them when necessary and simplify them 
when possible. In addition, we give a more detailed 
discussion of the conductance within the Time Depen- 
dent Current-Density Functional Theory (TDCDFT) for 
which we present a formally exact result. 

The theory is applied to molecular devices consisting 
of planar phenyl chains linked to gold wires via amine an- 
choring groups. We report theoretical values of the linear 
conductance for devices containing up to 4 phenyl rings, 
which are compared with the experimental data of Ref. 1 
Based on our analytic expression for Gc, we discuss and 
quantify the main factors influencing the charge trans- 
port in these devices. We recall that a similar study was 



recently completed for molecular devices involving amine 
linked alkyl chains-^^^ 



I. TRANSPORT: GENERAL CONSIDERATIONS 

We consider a charge transport experiment involving 
a device made of a molecular chain attached to metal- 
lic leads (see Fig. [ij. The system is driven by a small 
time oscillating electric field EiY\r,t), whose effects are 
treated in the linear response regime. The dc regime is 
obtained by letting the frequency of the oscillation go 
to zero. The existence of a steady state is implicitly as- 
sumed. 

Within the Time Dependent Current-Density Func- 
tional Theory (TDCDFT) an d line ar response regime, 
the current density is given byl^^'^ 



j(r,w) = 




where a^^ is the equilibrium Kohn-Sham conductivity 
tensor. A local density approximation expression for 
E|"(r, cu) is given in Ref. [T5l 

Ef = -V(/)r + -V(/.r° + (2) 
e e 

where is the linearized Hartree-exchange-correlation 
potential of the equilibrium DFT and Ej^° is the dynam- 
ical part of Ef , given by F,f"=--^VC, with ( the vis- 
coelastic stress tensor. In the linear regime: 

/dr'CT^^(r,r')Ef"(r') 

(3) 

= / dr' / dr"tT'^^ (r , r')^(r' , r")j (r") , 

where 

^„^(r,r')^ . (4) 
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jF(r, r') is understood as a matrix with elements 
Tafj(v,v') and matrix multiplication is understood be- 
tween a and P and between T and j. This leads to 



j(r) = / dr' / dv"[l - a"""- *T]-\vy) 
XCT'^"(r',r")V<?!)''^(r"), 



(5) 



which is an RPA type expression for the current density. 
Here 



(6) 



is the driving potential plus the adiabatic response of the 
electrons. 

The net current flowing through the device is given by 



1 = J^dS J dr' J dr"[l - a-^^ * .F]-i(r,r') 
xd-'^^(r',r")V(?!)''^(r"), 



(7) 



where S is an arbitrary transversal section. The potential 
drop that is measured by a voltmeter attached to the two 
ends of the device is given byP^ 



(8) 



and the linear conductance is defined as G = Note 
that the screening also contributes to the potential drop. 
Here, 4>i is the Hartree potential corresponding to the 
density perturbation ni. 



A. An exact expression for linear conductance 

We now show that A(/) can be pulled out of the com- 
plicated integrals in Eq. [7] For this, let us restrict the 
integral over dr" in Eq.[7]to a volume between two distant 
sections E_ and S+. We will later take these surfaces to 
infinity. Now, because 

Y,d^a^^r,r') = J] a;,a^^(r, r') = 0, (9) 

a p 

we have 

a^"(r',r")V"r^(r") = V"(T'^^(r', r")0°"(r"), (10) 

and we can transform the integral over r" in Eq. [t] in a 
surface integral: 



I = J^dSjdr'{J dS"- J dS") 
x[l - ct'^" * j£']-i(r,r')<T«^(r',r")?!)"''(r") 



We now chose the sections E± to be iso-surfaces of <^'"* 
in which case: 



(12) 



But once we pulled the potential out, the surface inte- 
grals no longer depend on the shape and position of the 
surfaces, a fact that follows from the property in Eq. |9] 
Therefore, we can deform S± into one single surface S" 
to obtain: 



x[l - CT-^s * j£']-i(r,r')o-''^(r',r"). 
At this point, let us write the explicit expression of < 



(13) 



r-(r) = 0r(r) + /dr'|^ 



-/dr' 



5n(r') 



0r 



|r-r'| 

ni(r'). 



(14) 



The density ni is localized near the junction, but its de- 
cay away from the junction can be rather slow. Due 
to the long range of the Coulomb kernel, the Hartree 
potential will take finite values at ±oo and will con- 
tribute to Acj)""^. The contribution from xc part was 
discussed in Refs. Tfl and fTSl, Here, it was pointed 
out that the common density functionals use semi-local 
exchange-correlation potentials in which case the kernel 
6v^^{r) /Sn{r') decays extremely fast with the separation 
|r — r'l and therefore the last integral in Eq. 
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: will vanish 

when r is taken at ±oo. The conclusion is that A^"** is 
in fact the potential drop measured by a voltmeter (see 
Eq.[8|. 

However, in same references it was also pointed out 
that functionals like those involving exact exchange lead 
to kernels Sv^'-'{r) /Sn{r') slowly decaying with |r — r'|, 
in which case the last integral in Eq. [14] will take finite 
values when r is taken at ±oo. In this case, we have to 
treat (pf^ as we treated E^^", in which case Tapir, r') has 
to be redefined as: 



+ 5a0n(r 



Sjpir') 



(15) 



This expression has to be computed at finite frequencies 
first, where one will use the relation rii = J^Vi, and then 
the limit uj^O has to be considered. Of course, in this 
case one can no longer use the local approximation of 

In either case, we arrive at the following formally exact 
expression of the linear conductance: 



X [1 - o-«s * J^rHr, r')d-«"(r', r") 



(16) 

Here, rj^ and denote the coordinates of two normal 
surfaces to the axis of the device. The position of these 
two surfaces can be taken arbitrarily. 



B. Linear conductance in adiabatic approximation 

The adiabatic approximation neglects the dynamical 
effects, which is equivalent to setting J- to zero. In this 



case, the expression for the Hnear conductance becomes: 

G - Jdr^l drl a-(r^, z; r^, z'). (17) 

We should point out that this expression also assumes a 
rapidly decaying kernel dv^^{r) / Sn(r') with the separa- 
tion |r — r'|. It is remarkable that, after the inclusion of 
electronic screening in A(f>, the expression for G remains 
formally identical to the one derived by Baranger and 
Stone^^for non-interacting electrons. 

In the rest of the paper, the conductance will be eval- 



uated using Eq. [17) a choice that is largely dictated by 
practical considerations. It amounts to implicitly assume 
that the dynamical effects are small, an hypothesis that 
we are unable to support with rigorous arguments. A 
previous numerical study found that dynamical effects 
play only a minor role,'221 but this study used only the lo- 
cal density approximation for E''*'" and considered small 
junctions, while here we focus on long molecular chains. 
We also leave the issue of non-locality of the xc poten- 
tial to future investigations. Although very interesting, 
these studies would be extremely challenging, particu- 
larly in the case of the large systems that are considered 
in this paper. 

For the nonlocal zz component of the conduc tivity ten- 
sor, we can work with the following expression! ^^ ! ^^ ! 



where 



AG,,(r,r') d', AG,,(r',r), (18) 



AG,,(r,r') = G,,+,5(r,r') 



5(r,r'), (19) 



and Ge is the Green's function Gf—{e-H) ^ of the Kohn- 
Sham Hamiltonian describing the equilibrium of the en- 
tire device: ^f=-|^V^-hKtf, Ktt = Vp, + v'''^''[n] with v^, 
being the ions' pseudo-potential and n the electron den- 
sity at equilibrium. As pointed out in Refs. f^TlandlT^ 
azz contains additional terms but they cancel out after 
the integrations in Eq.[T7|and therefore can be neglected. 



II. AN ANALYTIC EXPRESSION FOR THE 
TUNNELING CONDUCTANCE 

Consider a molecular device consisting of a long but 
finite periodic molecular chain (of unit cell b) attached 
to infinite metallic electrodes, like in Fig. [l] The ori- 
entation of the chain is along the z axis. We assume 
that a self-consistent Kohn-Sham calculation has been 
completed for the entire device. The effective potential 
of the entire molecular device V^tt is decomposed into a 
perfectly periodic piece, Vq, extending from -co to +oo, 
and a difference A^^V^ff — Vq- The periodic potential 
Vq is constructed by periodically repeating the effective 
potential between —b/2 and 5/2 at the middle of the 
chain. Our main assumption is that the potential dif- 
ference AV=V^i! — Vq rapidly decays to zero inside the 




L=Nb 




FIG. 1: Illustration of a typical device considered in this pa- 
per. The figure indicate the unit cell that is repeated period- 
ically to obtain the periodic potential Vb- It also defines the 
length L of the chain. 



periodic chain. In other words, we assume that, to a 
very high degree, the effective potential inside the chain 
is periodic. This assumption proved to be accurate for 
the systems we studied so far, including the phenyl chains 
studied in this paper. 

We regard the self-consistent Kohn-Sham Hamiltonian 
of the chain-|-leads as a periodic Hamiltonian, 

fP 

Ho = - — S/^ + Voir), Voir + bez) = Voir), (20) 
2m 

strongly perturbed by the potential AV^. The effective 
Hamiltonian of the entire system is then 



H = Ho + AVLir) + AVR{r), 



(21) 



where we divided AT^ into left and right parts. We as- 
sume that AVl.r decay fast to zero as we move away 
from the contacts. We demonstrate in the following that, 
based on an analytic expression for the Green's function 
corresponding to Hq, we can derive an analytic, non- 
perturbative expression for the Green's function of the 
entire device. This is somewhat complementary to the 
approach presented in Ref. [22l which views the devices 
as periodic leads perturbed by the junctions. 



A. Computing the Green's function for the 
periodic potential 

Let us first consider the Green's function Gg = (e — 
Ho)~^, with e outside the spectrum of Hq. To make 
the discussion more transparent, we recall that in 1 di- 
mension, the Green's function for a Hamiltonian of the 
form — j^^TT + y{x) can be conveniently written as: 



w^(^<,V'>) 

with a;< = min(a;,a;') and a;> = max(a;,x') and 
W{ip,(j))—ip(l)'-(j)ip'. Here, tp<{x) and ip>{x) are the so- 
lutions of the Schrodinger equation: 



2m dx^ 



V{x)\il){x) = eij}{x) 



(23) 
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FIG. 2: The generic shape of the Riemann surface of the 
bands for strictly one dimensional periodic systems (left) and 
for periodic chains in 3 dimensions (right). The figure also 
illustrates the contour F used in the main text. 



decaying to zero as a; — > — oo and a; ^ +oo, respectively. 
For a periodic system, the above expression reduces to: 



2m V'_fc(a;<)-0fc(x>) 



(24) 



where 4'k{x) is the Bloch function evaluated at the unique 
complex k with lm[/c]>0 for which the complex band en- 
ergy satisfies = e. To understand the simplicity of the 
above expressions, one should compare them with the 
formal expansion: 



E 



(l)„{x)*4>n{x') 
En - e 



(25) 



where {(/)„(x),e„} is the infinite sequence of eigenvectors 
and corresponding eigenvalues of the Hamiltonian. As 
opposed to Eq. |23] in Eq. |25] one has to compute a large 
number of wavefunctions and a truncation to n=N will 
generate 0{N) errors. The expression shown in Eq. 23 
is generally valid only in 1 dimension. We will show in 
the following, however, that for periodic Hamiltonians we 
can derive this expression using the Riemann structure 
of the bands. Since the molecular chains in 3 dimensions 
still exhibit a Riemann structure,'22l such derivation al- 
lows us to generalize Eq. 1 24 [ from strictly 1 dimension to 
molecular chains in 3 dimensions. 

We start now the derivation. From Ref. i24i it is known 
that the Bloch function ipxix) and the band energy e\ 
[X = exp(ikb)] can be defined on a Riemann surface that 
looks like in Fig. [2] This Riemann surface is made of 
a sequence of unit disks that are cut and then re-glued 
together as explained in Ref. [21] Different disks corre- 
spond to different bands and the physical, real k bands 
can be generated by evaluating e\ along the unit circles 
of each disks. Starting from the eigenvalue-eigenvector 
expansion, where we use the standard normalization of 
the Bloch functions: 



1 



fc/2 



6/2 



'4!n,-k{x)i^n,k{x)dx = 1, 



(26) 



we can write: 

G,{x,x') = 



^^VV-feMVVfeM^ (27) 



By using the Riemann structure, we can combine the sum 
over the band index and the integration over k into one 
single integral over a contour F defined on the Riemann 
surface of the bands (see Fig. l2l 



G,(x,x') 



dX ipi/x{x)ilJx{x') 



2iTbX 



(28) 



Now note that by changing the integration parameter 
from A to 1/A we interchange x and x' . Since A and 1/A 
run over the same path F we can write 



G,{x,x') 



dX ■'Pi/\{x<)^x{x^) 



2'KbX 



(29) 



We deform now the contour F towards the origin. No- 
tice that contour goes smoothly over the branch points 
since F has components on each pair of Riemann sur- 
faces connected by the branch points. Also, when the 
contour nears the origin, the integrand goes to zero be- 
cause 4>i/\{x<)4'\{x >) converges to A'"^""^ thanks to 
the correct ordering of x and x' . Thus, the only sin- 
gularity encountered during the deformation process is 
when e\ brushes over e and from the Residue Theorem 
we obtain: 



G,{x,x') = 



^i/\{x<)-ip\{x>) 
ibXdxex 



(30) 



If we go back to the k representation, the above expres- 
sion is the same as the one written in Eq. [24] and this 
ends our proof for the strictly one dimensional case. 

For periodic molecular chains in 3 dimensions, the Rie- 
mann surface of the bands was discussed in Ref. [531 and 
a typical shape is shown in Fig.|2]D. The difference is that 
now on each disks we can have more than two algebraic 
branch points and the equation = e has an infinite 
sequence cx^ of solutions. Starting from the expression 



Ge(r,r)=^ 



dX V^i/A(r<)-0A(r>) 
27rA ea ^ e 



(31) 



where F is the contour shown in Fig. [2Jd, and deforming 
the contour towards the origin and applying the Residue 
Theorem we obtain: 



G,ix,x')^J2 



V'i/A„(r<)V'A„(r>) 



iXadxex^ 

In the k representation, this expression becomes: 

V'-fe„(r<)V'fc„(r>) 



G°(r,r') = 5: 



(32) 



(33) 



where {fca} is the infinite sequence of wavenumbers with 
Im[A:]>0 such that e^^ = e and r</r> r/r' if z < z' 
and r</r> = r'/r otherwise. 
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B. Computing the Green's function for the entire 
device 



We can show in just a few steps why Eq. [33] is useful. 
Indeed, the Green's function for the entire device: = 
{H — e)^^ can be computed from the identity: 



G,(r,r') = GO(r,r') + /dr"/dr" 
xGO(r,r")T,(r",r"')GO(r"',r'), 
where the matrix is given by 

Te = + AVGeAV. 



(34) 



(35) 



Given that AV—AVl+AVr, we can naturally decom- 
pose the T matrix as 



(36) 



Now, the key observation is that, because of the localiza- 
tion properties of AVl and AVu, by taking r and r' near 
the middle of the chain we can tell what is the ordering 
between r and r" and between r'" and r' in the inte- 



gral of Eq. 34 Given this, the integrals can be formally 
executed and the result is: 



G,(r,r') = G«(r,r')+E J 



{r^^^k^ (r)V'fc, (r') + T^^^-k^ {v)i^-k, (r') 
+r,f Vfc.(r)V'-fe,(r') + Tnf V_fc„(r)^fe,(r')} 



(37) 



where 



na/3 



/dr/dr' i/._fcJr)r,(r,r')V-fc,(r') 



= JdvJdv' ^feJr)T,(r,r')Vfe,(r') 

^ LR 



Jdvjdr' V-fc„(r)r,,(r,r')^fe,(r') 



(38) 



T^,^ = / dr / dr' (r)r,,(r, r')^-;o, (r') 



This is the analytic expression of the Green's function we 
mentioned at the beginning. The "T" coefficients remain 
to be computed numerically, but at this point we have 
obtained the exact dependence of G^ on the coordinates r 
and r', which will allows us to compute the conductivity 



tensor. Eq. 33 is also essential for deriving the exact 
asymptotic form of the "T" coefficients in the limit of 
long chains.'^ 

Since Hq has no spectrum at ep, G^ behaves smoothly 
when e crosses the real line and consequently (see Eq. 19 1: 



AG,,(r,r')=E 



^ idk^kc'^&k^kr^ 
OL,f3 



[AT^^4,,^{r)^k,{r') + Ar,"'^V_,„(r)V^_,.,(r') 
-f AT,"/i^fe„ (r) V-fc, (r') + AT^U.^^ (r) V^fc, (r') } 
where AT stands for T^^^is — T^^^ig. 



(39) 



C. Computing the tunneling conductance 



Given our expression for the Green's function Eq. 39 



it 

is evident that the integrals in Eq. [T7|lead to generalized 
Wronskians between different Bloch functions. The gen- 
eralized Wronskian for two functions if) and (j) is defined 
as: 



M^(V',0)= / dr_L i/'(r_L,^;)9^(/)(r_L,2;). 



(40) 



We have the following remarkable propertyJTH valid at 
arbitrary energy e: 



(41) 



where {ka} is the sequence of wavenumbers correspond- 
ing to the energy e. Applying the above rules, we obtain 
the following expression for conductance: 



G 



Q,/3 



idkCkjOkek 



(42) 



The above expression is exact for insulating chains. It 
does not apply to metallic chains since we used the fact 
that Hq does not have spectrum at the Fermi level. The 
matrix elements of AT have simple and intuitive expres- 
sions: 

AT^'' = J dr J dr' 
x^A-fe„(r)AK(r)AGe,(r,r')AK(r') 

AT^^'^^/rfr/dr' Vfc„(r) 
X V'-fc, (r') AK.(r) AG,, (r, r') Ay,(r') Vfc, (r') 

(43) 

AT^J' ^ f dr J dr' 
X V'-fe. (r) AK(r) AG,,(r, r')AK,(r')^fe, (r') 

AT,f = /dr/dr' 
x^k^ (r) AyH(r) AG,,(r, r')AV,{r')^^k, (r') 

and they can all be expressed in terms of the spectral 
operator Pef. = ^i[G^+ - G^-]. The diagonal Pe,(r,r) of 
the spectral operator gives the local density of states. 

As discussed in [12] AT^f and AT,^^ coefficients are 
exponentially small compared to AT^^ and AT"^ , so in 
the asymptotic limit of long molecular chains, the tun- 
neling conductance is given by: 



2e2^ e'^'^e'^'^ 



J{ka+kp)L 



h ^ idkikjdkekn 



(44) 
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FIG. 3: Atomic configurations of the molecular devices. Diflterent rows show different view angles. 



with: 



and 



er = 2^/dr/drV_feJri,z+f)x 
AK(r)p,,(r,r')AK(r')^-fc,(rl,z'+f) 



e 

(48) 

This is the way we actually compute the coefficients in 
this work and details about how we choose 6 will be given 
later in the paper. An alternative way will be to compute 
the spectral operator directly from the Green's functions. 
This involves inverting the large matrices {H — ezL6)~^, 
which can be done iteratively and would not require sav- 
ing large amounts of data. 



(45) 



er = 2nJdvJdr'^kAr^,z-^)x 
AK(r)pe,(r,r')AK(r')^fe,(rl,z'-f). 



In the limit L oo, the Q coefficients become inde- 
pendent of L. Strictly speaking, the asymptotic form of 
G{L) is determined by the wavenumber k with minimum 
imaginary component. This is the case for the phenyl 
chains that we will investigate in the next Section, or for 
alkyl chains that were investigated in Ref. 13. However, 
for more complex molecular chains such as carbon nan- 
otubes ESl there may be many wavenumbers with similar 
imaginary parts, in which case we must consider more 
than one evanescent channel in Eq. 44 We point out 



that Eq. |44j tells how the evanescent tunnels interfere 
with each other during tunneling transport. 

It is important to observe that computing the contact 
conductance requires a converged density of states near 
the contacts, which can be obtained from a standard su- 
percell calculation that includes large enough electrodes. 
The spectral operator p^^ (r, r') can be computed in var- 
ious ways and each way can have its advantages and 
disadvantages. Provided one can store a large number 
of orbitals, a straightforward way consists in using the 
Kohn-Sham orbitals (A,: 



which leads to: 



lj: ^,.J^.^S^ ^:(r:)Mr% (46) 



e 



2(5 



J dr ^_fcJri,z+f)AK(r)0:(r) 



/dr0,(r)AK(r)^A-fc^(ri,2+f) 



(47) 



III. APPLICATION TO DEVICES INVOLVING 
PHENYL CHAINS 



In the following, we present an application to devices 
made of phenyl chains attached to gold electrodes via 
amine groups, like the those investigated in Ref. [H The 
complex band structure calculations of Ref. i26] reveal an 
evanescent channel with Im[/c] much smaller than that of 
the rest of the channels. Consequently, the tunneling con- 
ductance is determined by this evanescent channel and 
the expression for the tunneling conductance simplifies 
to: 



G - e^eHe- 



2ikL 



(49) 



This is to be compared to the classical expression 
G—Gce~^^. The tunneling coefficient (3 is related to 
k via (3 = 2Im[fc]&. The contact conductance Gc is given 
by the pre-exponential factor in Eq. 49 To be precise. 



let us write the simplified expression of theta coefficients 



2tt 



V'-fc(r)AK(r)p.,(r,r')AK(r')V'-fe(r'), 



(50) 
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FIG. 4: A plot of the lateral average of the effective potential 
of a 5 layer Au slab, corresponding to various values of pa- 
rameter Q. The vacuum region around the slab is much larger 
than what is visible in the picture. 



with r and r' measured from the left end of the chain. 
Similarly 

©H = m^fcy/rfr/dr'x ^^^^ 

with r and r' measured from the right end of the chain. 
We have included the derivatives idket into the 6 coef- 
ficients, and then we expressed this derivatives using the 
generalized Wronskian. This particular way of writing 
the Q coefficients is useful since the formulas become in- 
dependent of the normalization of the evanescent waves. 



A. Computational details 

We study three devices, containing 2, 3 and 4 phenyl 
rings, linked to gold electrodes via amine groups. These 
three devices will be referred to as (a), (b) and (c), re- 
spectively. The corresponding atomic configurations are 
shown in Fig. |3] Only the planar configuration for the 
phenyl chain will be considered. The geometry of the pla- 
nar phenyl chain was build from the structure of biphenyl 
molecule reported in Ref. HT] This reference reports an 
average C-C bond length of 1.40 A for the ring C atoms 
and a separation between the phenyl rings of 1.49 A. 
With this bond lengths, the unit cell size of the chain is 
4.315 A in z direction. The C-H bond length was fixed at 
1.10 A. The bond lengths reported in Ref. [57] are weakly 
dependent on the functional and basis set being used in 
the calculations. 

The bond angles for the N atoms of the linking groups 
were fixed in a tetrahedral configuration, except for the 
bond with the Au atom. The N-C and N-H bond lengths 
were fixed at 1.41 A and 1.04 A, respectively. The Au-N 
bond length was fixed at 2.40 A and the C-N-Au bond 
angle was fixed at 123°. Indicating by A, B, and C the 
stacking planes in the (111) direction for fee Au, the de- 
vices can be represented schematically by: 

CBACBA-Au-NH2-(C6H4)Ar-NH2-Au-CBACBA (52) 



Ideally, the left (right) Au ad-atom would occupy a lat- 
tice site of the C (A) stacking plane. Because of compu- 
tational constrains that require the chain to be oriented 
along the z direction and the surface of the electrode to 
be perpendicular to the z direction, this ideal configu- 
ration cannot be exactly satisfied, instead the ad-atoms 
are displaced towards the chain's plane by about 0.5 A. 
Since we are interested here mainly in illustrating the 
method, we did not investigate the issue of how geomet- 
rical factors such as the accurate position of the adatoms 
and distortions of the phenyl chain affect the calculated 
conductance. These issues are, however, very important 
for accurate quantitative comparisons with experiment. 

The lattice constant for the gold atoms in the leads 
was fixed at the experimental value (thus the stacking 
planes are spaced by 2.35 A). No surface reconstruction 
was considered. The system is periodically repeated in all 
three direction, but the calculations are restricted to the 
r point. In the z direction, the periodically repeated sys- 
tem has 12 layers of Au between two consecutive phenyl 
chains. For such electrode size, we expect the density of 
states near the contacts to be well converged. The lateral 
size of the supercell was chosen so that 20 Au atoms are 
contained in each layer. Thus, our computational super- 
cell contains 242 Au atoms. In total, there are 268, 278 
and 288 atoms for devices (a), (b) and (c), respectively. 

The equilibrium self-consistent Kohn-Sham calcula- 
tions were performed with a real space, pseudopotential 
code based on finite differences. The same code was used 
for the calculations reported in Ref. [T31 We adopted a 
5-point finite difference approximation for the kinetic en- 
ergy operator, and used a uniform rectangular space grid 
with a spacing of 0.3547 a.u., sufficient for a good conver- 
gence of the electronic structure. This grid is commen- 
surate with the unit cell of the periodic phenyl chain, 
which is the reference system in our transport calcula- 
tions. We adopted the Local Density Approximation 
(LDA) for exchange and correlation using the Perdew- 
Zunger (PZ)^^ interpolation of the numerical electron-gas 
data of Ceperley and Alder.[22l We used TrouUier-Martin 
norm-conserving pseudo-potential^^ for all the atomic 
species. The pseudopotentials for C and N atoms had 
distinct s and p components and we took the p pseudo- 
potential as the local reference. Purely local pseudopo- 
tentials were used for the H and Au atoms. In the latter 
case only the outermost s electrons were treated explic- 
itly. 

Since the current calculations include only the s elec- 
trons of Au, the calculated work function of the leads 
differs from the experimental value. To address this prob- 
lem, non- linear core corrections were proposed in Ref. 1311 
Even with these corrections, the work function of fee 
Au, as given by LDA calculations, takes values between 
6.28 and 6.70 eV, depending on the surface orientation 
l32l On the other hand, when the d Au electrons are 
treated explicitly, the LDA yield^^ workfunctions close 
to experiment.^ In our calculation, with the non-linear 
core corrections implemented as in Ref. 31, we find a 
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FIG. 5: The planar average (over xy) of the local density of states, shown as a density plot with energy (in cV) on the vertical 
axis and z coordinate (in a.u.) on the horizontal axis. 



work function of 6.6 eV for the Au electrodes, which 
should be compared to an average experimental value 
of 5.4 eV for the workfunction of the Au (111) surface. ■^'^ 
While this difference had insignificant consequences for 
the alkyl chains, ^'^ due to their large insulating gap and 
to the particularities of their complex band structure, for 
phenyl chains the consequences will be more severe due 
to their smaller insulating gap and to the parabolic shape 
of the complex band. More precisely, the Fermi level of 
the device will be located extremely close to the edge of 
the valence band of the insulating chain. 

Since our main purpose here is to demonstrate our 
methodology, we adopted a simple empirical approach to 
correct this shortcoming: we modified the local pseudo- 
potential of Au atoms by adding a local core correction 
of the form and{r) (Ry), where ^^(r) is the density of 
the frozen Au d electrons. The work function for the Au 
(111) surface becomes 5.4 eV if the constant a is fixed 
at 2.5 RyxBohr^. Fig. |4] shows plots of the effective po- 
tential of a 5 layer Au slab for increasing values of a. 
Here we can see a monotonic bending of the potential in 
the vacuum region, leading to a reduction of the work- 
function. We can also see a relatively large change inside 
the d cores, but these changes have minor effects on the 
occupied electron density since they occur well above the 
Fermi level. In addition, we do see a small change in the 
potential in between the planes. 



B. Electronic Structure 

The results of the electronic structure calculations are 
summarized in Fig. |5] which illustrates the local density 
of states for the three devices, averaged in the xy plane: 
pg^y{z, e)—J p^(x,y, z)dxdy. The plots give a color map 
of yOav(-z, e) in the plane of energy e and of position z. The 
figure was constructed from all Kohn-Sham orbitals used 
in the transport calculations, their number being equal to 
the number of the occupied orbitals plus additional 110 
un-occupied orbitals (without counting the spin). The 



Fermi level was fixed at zero and is indicated by the red 
line. In these plots, the conducting states of the leads and 
the band edges of the insulating chain are quite visible. 
The Fermi level, which is pinned by the continuum states 
of the leads, falls into the insulating gap of the phenyl 
chain. One sees that the conducting states of the leads 
decay rapidly to zero inside the phenyl chain, where the 
spectral gap becomes visible. The gap is clean all the 
way to the first gold atoms of the electrodes, showing no 
surface resonances. For energies inside the spectral gap 
of the chain, pay{z, e) does not show any special features 
near the contacts. The insulating band gap seen in Fig. [5] 
is larger for device (a) and is comparable for devices (b) 
and (c) . The tunneling transport is sensitive to the Fermi 
level alignment relative to the edges of the insulating gap. 
We point the reader to the Refs. '34*35' which give an 
extended discussion of the band alignment in molecular 
electronic devices and its effect on transport. 

Fig. |6] illustrates the local part of the AV, confirming 
the main assumption behind our formalism, namely that 
the potential inside the insulating chain is, to a very high 
degree, periodic and that is localized on the leads. 
Since we use norm conserving pseudo-potentials, the non- 
local part of AV is automatically localized on the leads. 

The band structure of the periodic potential Vq for 
device (c) is shown in Fig. |7] The real and complex 
structures are similar to those reported in Ref. ESI at 
least for energies below the vacuum. Above the vacuum, 
our calculation shows additional bands originating from 
scattering states, which are absent in the tight binding 
calculations of Ref. 



C. The conductance: Numerical results. 

We would like to comment first on the numerical ad- 
vantages brought in by our formalism. Due to the large 
supercells involved in this kind of calculations, very of- 
ten transport calculations for long molecular chains are 
carried out in a reduced basis set representation of the 
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FIG. 6: a) Atomic configuration of device (c). b) An iso- 
surface of Vgg-. c) Planar average of Vggf (with respect to the 
xy coordinates), d) An iso-surface of AV. e) Planar average 
of AV (with respect to the xy coordinates). The energy units 
are Ry. 



Hilbert space of the electron states. This can be problem- 
atic because the basis set functions are usually localized 
and it is not always clear how well are the scattering 
states represented by a small number of such functions. 
When considering experimental values for G that are be- 
tween 10"'^ and 10~^Go or even smaller, one can easily 
see that there is very little margin for errors. In our cal- 
culations, all the quantities involved in the formula for G 
are computed on the same grid used for the self-consistent 
calculation. Since the asymptotic expression of G is vir- 
tually exact for long chains, the analytic formula of Eq.[44| 
allows us to compute G without truncating our Hilbert 
space. 

We computed the transmission coefficient of our de- 
vices by evaluating Eq. 



49 at several energies e within 



the insulating gap and the results are reported in Fig. [8] 
as a function of e-ep. We should point out that the com- 
puted values become less accurate for energies closer to 
the band edges. The calculated transmission of the de- 
vice (a) looks different from the others, mainly because 
of its larger insulating gap. The linear conductance of 
the three devices, as derived from these calculations, are 
G = 1.5 X 10-3, 1.5 X 10-3, 4.3 x IO-^Gq, respectively. 
The /3 coefficient, computed as 26Im[fcj;-], is equal to 1.15 
for device (a) and 0.98 for the other two devices. It ap- 
pears that only the last two devices reached the asymp- 
totic tunneling regime. However, the situation is highly 
dependent on the position of the Fermi level. For ex- 
ample, /3 would be the same for the three devices if the 
Fermi level would move away from the valence band edge 
of the phenyl chain by 0.2 eV. Since the values of G 
are highly sensitive to the band alignment, we should 
be cautious when comparing the theoretical predictions 
with the experimental values. In any case, the predicted 
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FIG. 7: Real (right panel) and complex (left panel) band 
structures corresponding to the periodic potential Vo for de- 
vice (d). Only the complex band with smallest Im[fc] is shown. 
The Fermi level of the device was set to zero. The energy unit 
is eV and the unit for Im[fc] is 1/Bohr. 



G for device (a) is very close to the value measured in 
Ref. lU The predicted value of device (b) is 8.3 times 
larger than the experimental value reported in Ref. [T] 
No experimental value has been reported for the device 
(c). It is interesting to remark that a previous study'SSl 
on a device consisting of a single phenyl molecule linked 
to gold electrodes via amine groups predicted a theoret- 
ical G that is 7 times larger than the measured experi- 
mental value. The same reference pointed out that the 
calculated DFT conductance would become comparable 
to the experimental value if the Fermi level were located 
0.5 eV further away from the valence band. We also see 
from our data that a shift of ep by 0.5 eV would bring 
the theoretical prediction for both devices (a) and (b) in 
line with the experimental values. 

We now describe how we computed the conductance. 
The complex band structure corresponding to Vq varies 
slightly when different devices are considered. Overall, 
the band structure for Vo is similar to that reported in 
Ref. 26 for the infinite, isolated phenyl chains, suggest- 
ing that the main difference between Vq and the effec- 
tive potential of the infinite, isolated chain is a rigid 
shift. Given the particular complex band structure of 
the phenyl chains, the tunneling conductance is deter- 
mined by just one complex band, the one with the small- 
est Im[fc] . This complex band is shown in Fig. [Tjfor device 
(c). It was obtained by varying continuously Im[A;] from 
to its maximum value, while keeping Re[A:]= 0. For 
each complex value of k, the spectrum of the k depen- 
dent Hamiltonian: 

Hk = -{W-ike,f + Vo + e-'''^'~''^V^on-ioc{r,r'), (53) 

with periodic boundary conditions a,t z — ±&/2, was cal- 
culated and its eigenvalues ordered according to their real 
parts: Re[en;]<Re[e2fc]< • • • • We focus, in particular, on 
the 14th and 15th eigenvalues ei4k and €15^ (which take 
real values, see Fig. [t]) and their corresponding evanes- 
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FIG. 8: Plots of the transmission as function of energy; green, 
blue and red colors are used for devices (a), (b) and (c), re- 
spectively. 



cent Bloch functions ipiik and V'lSfc- When lni[fc]=0, 614^ 
and eisfc coincide, respectively, with the top of the va- 
lence bands and with the bottom of the conduction bands 
of Vq. By increasing Im[fc], the two eigenvalues move 
towards each other until they become degenerate when 
k reaches the branch point at Im[A:]=0.08 Bohr^^. At 



different values of Im[fc], we evaluated Eq. 44 for both 
e=ei4fe and e=ei5fc, using the corresponding evanescent 
Bloch functions ipiik and V'lSfc to compute the & coeffi- 
cients via formulas |50] and [51] The spectral kernel was 
computed directly from the Kohn-Sham orbitals of the 
full device as previously explained. The coefficient S was 
fixed at 0.01 Ry. This value is about an order of magni- 
tude larger than the average energy level spacing of the 
Kohn-Sham orbitals near the Fermi energy. 



D. 



Insight into the transport properties of phenyl 
chains 



The analytic result of Eqs. 44] 50 and 51 allows us to 
point several key aspects of the tunneling transport of 
our devices. Since the formulas involve overlap integrals, 
the new insight is obtained by looking at each physical 
quantity entering in the expressions of the Q coeffien- 
cients. 

A plot of the local density of states (i.e. the diagonal 
part of the spectral operator) was already given in Fig. [5] 
and a plot of \AV\ was given in Fig. [6] Fig. [9] shows 
a plot of the evanescent Bloch solutions of the periodic 
Hamiltonian with potential Vq for device (c), evaluated 
at the Fermi level. These functions are a property of 
the periodic Hamiltonian only, but their spatial decay 
is fixed by the P coefficient, which depends on the level 
alignment as discussed earlier. The contact conductance 
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FIG. 9: Different angle views of an iso-surface plot (blue 
surface) corresponding to 5% of the maximum value of the 
evanescent Bloch function |i/)_fe(r)| at the Fermi energy for 
device (c). For illustrative convenience, the evanescent Bloch 
function |?/>_fc(r)| was truncated at the right end. For refer- 
ence, we also show an iso-surface plot of the effective potential 
(in red). The iso-surface plots for |i/)fe(r)| can be obtained by 
mirror symmetry relative to the center of the device. 



Gc depends on the overlap of these evanescent functions 
with other physical quantities, and a plot like the one 
in Fig. [9] allows us to assess quantitatively the contact 
region that is relevant to tunneling transport. 

A main factor in our transport calculation is the over- 
lap between the evanescent Bloch function ipz^ki'*^) and 



AVl/ri the quantity 



^-L/R^r) = VTfc(r)AK/R(r), 



(54) 



which is exponentially localized at the left/right contacts. 
As a consequence the spectral operator in Eq. [46] is only 
needed in a region near the contacts. A plot of ^'L/n (r) 
for device (c) is shown in Fig. 10 This plot allows us 



to understand how the different Au layers contribute to 
the contact conductance Gc-^^ From the data we extract 
that the contact Au atom and the next two gold layers 
contribute to ^'L/R(r) by about 75%, while the remaining 
25% comes from the remaining layers. This information 
tells us that the conductance of our devices is primarily 
determined by the first three layers of Au atoms, an in- 
formation that could be useful when designing molecular 
circuits based on phenyls. The spatial spread of ^'L/R (r) 



along the device seen in Fig. 10 is more extended than 
the one found for alkyl based devices. This implies that 
the conductance of the present devices is more sensitive 
to the geometrical and chemical configuration of the con- 
tact, or to the orientation of the molecule relative to the 
molecular wires. 



IV. CONCLUSIONS 

In conclusion, we presented an extended discussion of a 
previously proposed theoretical approach for off-resonant 
tunneling transport. We added details where necessary 
and we greatly simplified the derivation of the asymptotic 
expression for the tunneling conductance. In addition. 
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FIG. 10: An iso-surface plot of |^L/R(r)j, corresponding to 
2% of the maximum value of |'I'L/R(r)| and the planar aver- 
age of |*I'L/R(r)j (with respect to the xy coordinates). The 
horizontal axes of the graphs are aligned. For convenience, 
we also included the atomic configuration of the device. 



we reported a formally exact expression for the linear 
conductance derived within the TDCDFT. 

The application to devices involving phenyl chains re- 
vealed several interesting facts. Based on our data, it ap- 
pears that only devices (b) and (c) reached the tunneling 
regime. Experimentally, one observes that the conduc- 
tance of devices containing one, two and three phenyls 
obey the tunneling conductance formula G=Gce~^^ 



with same Gc and /3. This can be a coincidence or it can 
be the real fact. Only measurements on longer phenyl 
chains can clarify the point. 

We found that the transport calculations are extremely 
sensitive to band alignment. This is prompted in the first 
place by the relatively small insulating gap of the phenyl 
chain but also by the fact that LDA places the Fermi level 
close to the edge of the valence band of the phenyl chain. 
For this reason, kp is located in the rapidly varying region 
of the complex band and small variations in lead to 
large variations in conductance. 

The analytic expression for the tunneling conductance 
allowed us to probe several transport characteristics of 
the devices. We showed that the contact conductance 
is exponentially localized near the contacts and we were 
able to describe quantitatively this localization. Since 
the evanescent conducting channels decay slower than 
for the case of alkyl devices, the contact conductance 
is less localized and the tunneling characteristics of the 
phenyl based devices are predicted to be more sensitive 
to the particularities of the electrodes when compared to 
devices involving alkyl chains. 
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